Calculation of fission observables through event-by-event simulation 
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The increased interest in more exclusive fission observables has demanded more detailed models. 
We present here a new computational model, FREYA, that aims to meet this need by producing 
large samples of complete fission events from which any observable of interest can then be extracted 
consistently, including arbitrary correlations. The various model assumptions are described and the 
potential utility of the model is illustrated by means of several novel correlation observables. 
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I. INTRODUCTION 

Nuclear fission presents an interesting and challenging 
physics problem which is still, about seventy years af- 
ter its discovery, relatively poorly understood. Although 
much of the key physics involved is understood quali- 
tatively, a quantitative description is still not in sight, 
despite vigorous efforts by many researchers. 

Because of its inherent complexity, fission provides an 
important testing ground for both static and dynamical 
nuclear theories. Furthermore, fission is also important 
to society at large because of its many practical applica- 
tions, including energy production and counterprolifera- 
tion, topics of current urgency. 

Whereas the more traditional treatments of fission (see 
Ref. [l[ and references therein) have sought to describe 
only fairly integral fission properties, such as the average 
energy release and the average differential neutron yield, 
many modern applications require more exclusive quan- 
tities, such as fluctuations in certain observables (e.g. the 
neutron multiplicity) and correlations between between 
different observables (e.g. neutrons and photons). There 
is thus a need for developing models that include the 
treatment of fluctuations and correlations. 

A potentially powerful approach towards meeting this 
challenge is to develop simulation models that can gener- 
ate samples of complete fission events, since a subsequent 
event-by-event analysis could then provide any specific 
correlation observable of interest. Furthermore, due to 
the more detailed quantities that they can address, such 
models can provide valauble guidance to experimental- 
ists with regard to which observables are most crucial for 
further progress in the understanding of fission. 

Relatively recently, Lemaire et al. *2j presented a 
Monte-Carlo simulation of the statistical decay of fission 
fragments from spontaneous fission of 252 Cf and ther- 
mal fission of 235 U by sequential neutron emission. That 
work demonstrated how fission event simulations, in con- 
junction with experimental data on fission neutrons and 
physics models of fission and neutron emission, can be 
used to predict the neutron spectrum and to validate 
and improve the underlying physics models. 

We have developed a conceptually similar calculational 



framework within which large samples of complete fission 
events can be generated, starting from a fissionable nu- 
cleus at a specified excitation energy. The associated 
computational code is denoted FREYA (Fission Reaction 
Event Yield Algorithm) . We present here the model in its 
most basic form which, though quite simplistic in many 
regards, is already capable of producing interesting re- 
sults, as we shall illustrate. Furthermore, FREYA was em- 
ployed in a recent study of sequential neutron emission 
following neutron-induced fission of 240 Pu [3|. 

In its present early form, FREYA ignores the possibility 
of neutron emission from the nucleus prior to its fission 
(n th chance fission), and its applications are therefore 
limited to lower energies, such as thermal fission. 

In Sects. |n] and |TTT] we describe above how a single 
fission event is being simulated in the pilot version of 
FREYA. By repeating the procedure a large number of 
times, we may generate an entire sample of final fission 
events, each one consisting of two (slightly excited) resid- 
ual product nuclei and the various emitted neutron and 
photons, each one with its associated momentum. In the 
development of the numerical code, special care has been 
taken to design the various algorithms for fast execution. 
As a result, FREYA runs fairly fast, thus making it practi- 
cal to generate sufficiently large event samples to permit 
detailed correlation analyzes. In Sect. |IV] we discuss a 
number of illustrative results. 



II. FISSION 

When the possibility of pre-fission radiation is ignored, 
the first physics issues concern how the mass and charge 
of the initial compound nucleus is partitioned among the 
two fission fragments and how the available energy is di- 
vided between the excitation of the two fragments and 
their relative kinetic energy. 



A. Fission-fragment mass and charge distributions 

In our current understanding of the fission process, the 
evolution from the initial compound nucleus to two dis- 
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tinct fission fragments occurs gradually as a result of 
a disipative multidimensional evolution of the nuclear 
shape. However, since no quantitatively reliable theory 
has yet been developed for this process, we employ empir- 
ical evidence as a basis for selecting the mass and charge 
partition. Thus, the mass and charge partition of the fis- 
sioning nucleus a "Zq is determined by first selecting the 
mass partition from a specified probability distribution 
P(Af) and subsequently selecting the charge partition 
from the associated conditional probabibilty distribution 
P Af {Z f ). 

In a given event, the mass number Af of one of the 
fission fragments is selected randomly from a probability 
density P{Af) for which we employ five-gaussian fits to 
the product mass number distribution Q shifted upwards 
in mass to ensure a symmetric distribution of the primary 
fragments, 

m=+2 

P(A f ) = Yl ^mGm(Af) , (1) 
m=-2 

where each of the five Gaussians has the form 

Gm{Af) = (27T^p e -(A f -A f -D m f/2,l (2) 

Contrary to Ref. [ij], we are interested in the primary (i.e. 
pre-evaporation) fragment distribution rather than the 
final (post-evaporation) product distribution and there- 
fore use Af = \Aq. The fitted values of the normal- 
izations Af m — Af-m', the displacements D rn = D- m ; 
and the dispersions a m = <J- m have some dependence on 
the excitation energy Eg. Since J2a P{A) — 1 we have 
Afo + 2Nx + 2Af 2 = 1. 

It should be noted that the normalizations M m are 
not quite correct since the sum over the integer fragment 
mass numbers A f does not yield the exact integral of the 
Gaussian and the range of the fragment mass numbers 
Af is finite. Because neither of these inaccuracies plays 
a noticeable role, we shall ignore them in the present 
preliminary treatment. We also note that merely back- 
shifting the average Af but not reducing the widths a m 
(to take account of the smearing due to the neutron evap- 
oration) will lead to a product mass distribution that is 
a bit too wide (since the smearing effect of the neutron 
evaporation will, in effect, be taken into account twice). 
However, this effect is rather small and is ignored in the 
present treatment. 

For the subsequent selection of the fragment charge 
number Zf, we follow Ref. Q and employ a normal dis- 
tribution, 

P Af (Z f ) oc e -{z f -Zf)/2°l 7 (3) 

with the condition that \Zf — Z f\ < haz- The centroid is 
determined by demand that the fragments have the same 
charge-to-masss ratio as the fissioning nucleus, on aver- 
age, Z f — AfZf)/A . We use the values of the dispersion 
(Tz measured by Reisdorf et al. 5J, 0.40 for 236 U(n, /) 



and 0.50 for 239 Pu(n, /). [There appears to be an er- 
ror (presumably typographical) in the expression (2) for 
P(Z) in Ref. Q: the pre-exponential factor should be a 
square root in order for P(Z) to be normalized to unity.] 



B. Scission energetics 

We obtain the fission energetics by assuming that the 
two fission fragments lose contact at a certain scission 
configuration which we take to be two coaxial spheroidal 
prefragments with a specified tip separation d. For the 
time being, we ignore the nuclear proximity attraction 
between the two prefragments as well as any possible 
relative motion at the time of scission. These two effects, 
which counteract one another, are relatively small but 
should ultimately be considered. 

We introduce some degree of distortion of the pre- 
fragments relative to their ground-state shapes, due to 
their mutual Coulomb repulsion. This is done primarily 
in order to ensure that the resulting fragment excita- 
tions (and hence the neutron multiplicities) roughly re- 
semble those observed. Thus, generally, the deformation 
of the fragment at scission, e sc , is larger than that of 
the ground state, e gs . The associated distortion energy 
is calculated by using the small-deformation approxima- 
tion [i], SV = ^[E" - \E° c ]{el c - s 2 gs ), which suffices 
at this early stage of the development. (Here we use 
the macroscopic expressions for the surface energy Eg 
and the Coulomb energy Eq for the spherical shape, as 
described in App.[X]) The distortion moves the prefrag- 
ment centers apart, for any fixed tip separation d, and 
thus lowers the mutual Coulomb repulsion Vy. 

It follows that there are two contributions to total ex- 
citation of each prefragment, 

E* = SV + Q, , (4) 

namely the distortion energy SV% and the statistical ex- 
citaiton (heat) Qi. 

The Coulomb repulsion between the two deformed pre- 
fragments is calculated by means of the formula derived 
in Ref. [7[ for two coaxial, uniformly charged spheroids, 

V§ = e 2 *\, F{ Xi , Xj ). (5) 
J Ci + Cj + d 

The factor F is unity for two spheres and larger if one or 
both fragments are prolate. It depends on the dimension- 
less deformation measures Xi given by xf = (cf — bf)/ Rf , 
where Ci = Ri[l + — |e] 2 / 3 is the major axis and 

hi = Ri[l — + j^] 1 ^ 3 is the minor axis, while Ri is 

the average radius of the fragment. 

Once the fragments have lost contact, they are acceler- 
ated by their mutual Coulomb repulsion and their shapes 
relax to their equilibrium forms. The scission distortion 
energies are converted into additional statistical excita- 
tions of the respective fragments. We assume that these 
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processes have been completed before the de-excitation 
processes begin. 

With the (significant) simplifications described above, 
we have the following simple energy relations for any par- 
ticular fission channel, A °Z — ► Al Z l + Ah Z h , 

M* = Mq 8 + Eq = Mf + E* L + Mf + E* H + vEh 
= M* L + M* H + K LH . (6) 

Here Mf is the ground-state mass of the nucleus *Zi, 
i = 0,L, H, and E* is its excitation, so M* = Mf + E* 
is its total mass. [The ground-state masses are taken 
from the compilation by Audi et al.Jw , supplemented by 
calculated masses by Moller et al. Q where no data are 
available.] Furthermore, Vf H is the Coulomb repulsion 
between the two light and heavy fragments at scission. 
This energy is, by fiat, fully converted into relative kinetic 
energy of the two receding fission fragments, K\i- Thus, 
in addition to ignoring any possible post-scission dissipa- 
tion, we also disregard any angular-momentum effects. 
While these effect arc expected to be small, it might be 
of interest to include them at a later time. The Q-value 
associated with the particular fission channel is given by 

Qo^lh = Mf + E*-Mf-M^ = K LH + E* L + E* H . 

(7) 

C. Thermal fluctuations 

Once the scission configuration is known, its average 
total internal (statistical) excitation energy, Q, can be 
readily obtained, 

Q = Ql + Qh = Mf + E* - M S L C - MH ~ VEh , (8) 

where Mf c = Mf + SVi is the mass of the distorted pre- 
fragment of the scission configuration. We assume that 
this internal energy Q is partitioned statistically between 
the two prefragments, as would be the case when the two 
are in mutual thermal equilibrium. Thus, on the aver- 
age, the total excitation energy is divided in proportion 
to the respective heat capacities. These in turn are char- 
acterized by the Fermi-gas level-density parameters ai 
which are approximately proportional to the fragment 
masses Af, we use the values calculated in Ref. [l(| (see 
App. IB) . [We note that those calculations were made for 
nuclei in their ground-state shapes, whereas the scission 
prefragments are distorted and may thus have different 
effective level-density parameters.] The mean excitation 
in a nucleus is assumed to be Q i = a%Tf. so the heat 
capacity is dQi/dTi = 2a{Ei oc a,. Since the two pre- 
fragments in the scission configuration have a common 
temperature, T LH = [Q/{a L + a H )] 1/2 = [Q 4 /ai] 1/2 , we 
use Q i = a t Tl H . 

The fluctuations in the statistical excitation Qi are 
given by the associated thermal variances, of — 2Q(Elh- 
The fluctuations 5Qi arc therefore sampled from normal 
distributions with variances of . The prefragment excita- 
tions in a given event are then Qi = Qi + &Qi- 



As a result of the fluctuations in the statistical exci- 
tation energies of the individual prefragments, Qi, the 
combined statistical excitation energy, Q = Ql + Qh, 
will also fluctuate. This fluctuation in turn implies a 
compensating fluctuation in the total fragment kinetic 
energy, so that Klh — Klh + SKlh where 

Klh = V£ H , SK LH = SQ L - SQ H . (9) 

We note that the resulting thermal distribution of heat 
in each prefragment is approximately gaussian, 

Pi(Qi) « (27r ?)-ie-W«-^) a /a^ . ( i ) 

Consequently, the distribution of the combined amount 
of heat in both fragments, Q = Ql + Qh, is also approx- 
imately gaussian and the associated variance is the sum 
of the individual variances, Oq = a\ + cr H . Energy con- 
servation implies that the distribution of the total kinetic 
energy Klh is a gaussian with the same width, ok = ctq, 
as was assumed in Ref. [2] . 

It is physically reasonable that the partioning of the 
total energy between kinetic energy and internal excita- 
tion fluctuates because the evolution of the fissioning sys- 
tem from saddle to scission is a dissipative process. The 
associated conversion of the collective energy to heat is 
the result of many elementary stochastic processes. The 
fluctuation-dissipation theorem then relates the average 
energy loss (the dissipation) to the associated fluctua- 
tion. Energy conservation demands that the fluctuations 
in the kinetic energy are exactly the opposite of those 
in the internal excitation. These, in turn, are given by 
the above thermal expressions insofar as statistical equi- 
librium is maintained during the shape evolution from 
saddle to scission. [We ignore the possibility that the 
scission configuration itself might also fluctuate from one 
event to another for a given fission channel.] 

Once the relative kinetic energy Klh has been ob- 
tained as described above, the magnitude of the rela- 
tive momentum, plh, of the fully accelerated fragments 
is then determined. Since the kinetic energy is rela- 
tively small (Klh « 200 MeV, while M * > 200 GeV), 
we may safely assume that Klh *C M* and use non- 
relativistic kinematics, p\ H = 2^lhKlh, where the re- 
duced fragment mass is hlh = MlM H /{Ml+M H ) with 
M* = Mf c + Q t = Mf + 5Vi + Qi being the total mass 
of the excited prefragment. Ignoring any angular mo- 
mentum effects, we select the fission direction V ran- 
domly. The fragment momenta are then Pl = PlhV 
and P h = —PlhV , in the frame of the fissioning nu- 
cleus. 



III. POST-FISSION RADIATION 

As mentioned above, we assume that the two excited 
fragments do not begin to de-excite until after they have 
been fully accelerated by their mutual Coulomb replusion 
and their shapes have reverted to their equilibrium form, 
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which we take to be those of their ground states. [In prin- 
ciple, the equilibrium shape of a nucleus depends on its 
excitation since both shell effects and surface tension are 
temperature dependent, but we have ignored this rela- 
tively minor complication at this time.] Furthermore, we 
ignore the possibility of charged-particle emission from 
the fission fragments. 

Each of the fully relaxed and accelerated fission frag- 
ments typically emits one or more neutrons as well as 
a (larger) number of photons. We assume that neutron 
evaporation has been completed (i.e. no further neutron 
emission is energetically possible) before photon emission 
sets in. This simplifying assumption obviates the need for 
knowing the ratio of the widths, T 7 (E*)/T n (E*). 

A. Statistical evaporation of neutrons 

We treat post-fission neutron radiation by iterating a 
simple treatment of a single neutron evaporation, until 
no further neutron emission is energetically possible. 

Statistical neutron evaporation is but one example 
of a general two-body decay. In the present case, the 
initial body is an excited nucleus with a total mass 
equal to its ground-state mass plus its excitation energy, 
M* = Mf a + E* . The Q- value for neutron emission is 
then Q n — M* — Mf - m n , where Mf is the ground- 
state mass of the daughter nucleus and m n is the mass 
of the (unexcitable) ejectile (the neutron). The Q- value 
equals the maximum possible excitation energy of the 
daughter nucleus, which is achieved for vanishing final 
relative kinetic energy, Q n = i?™ ax , which would be ob- 
tained if the emitted neutron had no kinetic energy. It 
is related to the associated maximum daughter temper- 
ature T} nax by a f {Tp x ) 2 = £} nax , where a/ is the level 
density parameter of the daughter nucleus (see App. |B|) . 

1. Spectral profile 

Once the Q-value is known, it is straightforward to 
sample the kinetic energy of an evaporated neutron, as- 
suming that it is isotropic in the rest frame of the emit- 
ting nucleus. We first note that the kinetic energy of 
the neutron has the form e n = p 2 /2m n while v n oc ^Je^ 
(non-relativistically) so that d 3 p n oc ^/e^de n for isotropic 
emission. The differential distribution is then [ll|, [12] 

~3 — d 3 p n cx y^e~ ea/T T^ y/e^de n dn 
a p n 

= e n e- e " /T r x d en( «7 , (11) 

in the rest frame of the emitting nucleus. The form 
y / e^'exp(— e n /T) can be understood as the product of the 
thermal occupancy of the neutron, oc exp(— e n /T), and 
its normal speed v n oc y/e which introduces a bias in 
favor of those neutrons that are moving perpendicular to 
the nuclear surface. 



The kinetic energy of the evaporated neutron, e u , is 
sampled by means of a specific fast algorithm that is de- 
scribed in App. [C] We note that the form of the energy 
spectrum implies that the evaporated neutron has a mean 
(relative) kinetic energy of (e n ) = 2T™ ax and an associ- 
ated variance of 2(T™ ax ) 2 . These expressions apply to 
the particular stage of the evaporation chain. Generally, 
the first neutron evaporated from the fragment will tend 
to have a higher energy than the second one, and so on. 

2. Kinematics 

Although relativistic effects are very small, we wish to 
take them into account in order to ensure exact conser- 
vation of energy and momentum, which is convenient for 
code verification purposes. We therefore take the above 
sample value e to represent the total kinetic energy in 
the rest frame of the mother nucleus, i.e. it is the kinetic 
energy of the emitted neutron plus the recoil energy of 
the residual daughter nucleus. The excitation energy in 
the daughter nucleus is then given by 

Ef = Qn - £n • (12) 

Since relativistic mass of the daughter nucleus is MJ = 
M| s + Ef , it is possible to calculate the momenta of the 
emitted neutron and the excited daughter as follows. 

Generally, if a particle of mass M decays into two par- 
ticles of masses mi and ma, those two particles are emit- 
ted back-to-back in the rest frame of the initial particle, 
with their momenta having equal magnitudes. Denoting 
this common momentum magnitude by p, application of 
elementary energy conservation yields 

M = E 1+ E 2 = [ml+p 2 ] 1 ' 2 + [ml+p 2 ] 1 ' 2 , (13) 

from which the magnitude p can be readily obtained, 

4M 2 p 2 = [M 2 -{m 1 +m 2 ) 2 ][M 2 -(m 1 ~m 2 ) 2 ] . (14) 

The individual energies, Ei — [p 2 + mf] 1 / 2 , may then 
be obtained subsequently. We employ the above formula 
with M = M*, mi = m n and m 2 = Mf = M? + Q n — e n . 

Assuming that the emission is isotropic (which follows 
from the neglect of angular-momentum effects), we may 
readily sample the direction of relative motion (1?, ip). 
The momentum of the ejectile is then 

p n = (p cosip sin^,p sirup sim?,p cos$), (15) 

while the recoil momentum of the residue is the opposite, 
Pf = —p n . These momenta are in the two-body CM 
frame, the frame of the mother nucleus, which would 
generally be moving. We therefore need to boost these 
momenta to the overall reference frame (see App. [Pi). 

The emission procedure described above may be re- 
peated until no further neutron emission is energetically 
possible. That happens when E*j < S a , where is 
the neutron separation energy for the daughter nucleus, 
S n = M( A Z) - M( A ~ 1 Z) — m n . 
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B. Statistical emission of photons 

Although, at this initial stage, our main focus is on 
neutron evaporation, we wish to also include an approx- 
imate treatment of photon emission. For this purpose 
we disregard nuclear structure effects and treat the post- 
evaporation photon cascade in a manner that is similar 
to the neutron emission described above. Clearly, this 
part can be refined by taking account of the specific level 
structure in the fission fragments. Because the photon is 
massless, we introduce an energy cutoff (see below). 

Furthermore, the vanishing photon mass causes it to be 
ultrarelativistic with p 7 c = e 7 and v 1 = c. Consequently, 

oc e^e-^fde^dn , (16) 

<2 P 7 

as was also used in Ref. [Hj]. For the first photon to 
be emitted, Tf is the temperature in the nucleus right 
after the last neutron was evaporated, dfTj — Ej, and 
generally it is the temperature before the next photon is 
emitted. 

The photon energy e 7 is sampled by a fast algorithm 
(see App.[C| and the nuclear excitation energy is reduced 
correspondingly, (E%)' = Ef — e n . The spectral shape 
(TTH|) yields an average photon energy of (e 7 ) = 37/ and an 
associated variance of 3T| , for a fixed value of Tf. Since, 
in principle, the continuous form of the spectrum leads to 
an infinite number of ever softer photons, we keep track 
of only those with an energy above a specified threshold, 
£ mm _ 200 keV. For photons above that threshold, the 
emission direction is sampled uniformly over 47r and a 
Lorentz boost is performed to express the emitted photon 
and the nuclear residue in the overall reference frame. 

This procedure is iterated until the nuclear excitation 
energy falls below the specified minimum value e™ m . 

IV. ILLUSTRATIVE RESULTS 

Here, we wish to illustrate the utility of single-event 
models like FREYA by presenting a number of correla- 
tion observables that could not be addressed with ear- 
lier codes which have tended to focus on more inclusive 
quantities. Obviously, the present preliminary version of 
FREYA involves a number of simplifying approximations 
and, consequently, the results cannot be expected to be 
numerically accurate. Certainly, for the most common 
observables, such as average multiplicities and spectra, 
the most reliable results can undoubtedly be obtained 
from the well-tuned codes that have long been available. 
We expect that event simulation codes will, in due course, 
achieve a similar level of accuracy. Meanwhile, they may 
serve as useful supplements with which is will be possible 
to address more detailed observables on an approximate 
level. 

While the main purpose here is to illustrate the kind 
of novel information that can be accessed with FREYA, we 
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FIG. 1: The mean total kinetic energy, K to t, of the two fission 
fragments, and the associated dispersion (bars), as a function 
of the mass number of the heavy fragment, Ah, for 0.53 MeV n 
on 235 U (bottom) and 239 Pu (top). The data from Nishio 0] 
(with a few representative dispersions), Tsuchiya [l5| ). and 
Wagemans are shown. The dispersions reflect the width 
of the kinetic energy distribution and are not (experimental 
or theoretical) uncertainties. [Color online.] 
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FIG. 2: The mean kinetic energy of a single fragment and the 
associated dispersion (bars) as a function of its mass number 
Af, for 0.53 MeV n on 235 U (bottom) and 239 Pu (top). Pu 
data from Nishio [l4| and Tsuchiya [l5[ are shown. [Color 
online.] 



wish to first show a number of more familiar observables. 
Throughout we consider fission induced by thermal neu- 
trons on 235 U and 239 Pu. Fission induced by higher- 
energy neutrons is not considered, since the possibility 
of pre-fission neutron emission (and the associated n th 
chance fission) has not yet been included. 
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A. Fission fragments 

The most basic observable is perhaps the product mass 
distribution P(A p ) which, by design, matches the fits to 
the observed data and thus need not be displayed. 

We therefore start by considering the kinetic energies 
of the fission fragments. Figure [T] shows the combined 
kinetic energy of both fragments, -Ktot, as a function of 
the mass number of the heavy fragment, Ah, while Fig. [2] 
shows the kinetic energy of a single fragment as a function 
of its mass number Af. These results exhibit the general 
observed features, though the detailed behavior is not yet 
expected to be accurate. 

The figures show the mean values of the kinetic ener- 
gies as well as the associated dispersions. A quick com- 
parison of the two figures suggests that the variance of 
the total kinetic energy is generally larger than the sum 
of the variances of the individual kinetic energies. This 
striking feature is an elementary consequence of momen- 
tum conservation. Since the two fragments emerge with 
opposite momenta, the fluctuations in their kinetic ener- 
gies are closely correlated. As a result, the sum of the 
variances of the two individual fragment energies, K L 
and Kh, is significantly smaller than the variance in the 
combined fragment energy Klh = Kl + Kh, namely 
a 2 {K L ) + a 2 (K H ) = [1 - 2A L A H /A 2 ] a 2 (K LH ). In par- 
ticular, for a symmetric split, Al — Ah, we have a(Ki) = 
\a(K LH ) hence a{K LH ) 2 = 2[(a(K L ) 2 + a(K H ) 2 }. 

While the total excitation of the emerging fragments 
is related to their total kinetic energy by energy con- 
servation, its partition is less straightforward, depending 
both on the relative heat capacities {i.e. level densities) 
and the scission fluctuations. Figure [3] shows the mean 
fragment excitation E t together with the associated dis- 
persion, as a function of the fragment mass number Af. 
In the present model, the division of the available energy 
between kinetic and excitation is sensitive to the degree of 
distortion of the scission pre-fragments, a property that 
in turn depends on the shell structure of the specific nu- 
clides. 



B. Neutron multiplicities 

The fission fragment excitation energies E*(Af) (see 
Fig. [3]) largely determine the multiplicities of evaporated 
neutrons v(Af). This correspondance is clearly seen in 
Fig. 0] which shows the mean neutron multiplicity V(A f ) 
and the associated dispersion a v (Af). We note that the 
observed sawtooth shape is roughly reproduced, though 
the detailed behavior is not completely satisfactory. 

The overall neutron multiplicity distribution P{y) is 
shown in Fig.[5l This figure also shows the separate mul- 
tiplicity distributions P{vl) and P{vh) for the number 
neutrons v L and vh that were emitted by the light or 
the heavy fragment, respectively, a quantity that is dif- 
ficult to obtain experimentally. The associated average 
multiplicities are shown in Table [I] (77.L = (vl), etc.), 
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FIG. 3: The mean excitation energy E (curves) of a fission 
fragment and the associated dispersion (bars), as a function 
of the mass number Af for 0.53 MeVn on 235 U (bottom) and 
239 Pu (top). [Color online.] 
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FIG. 4: The mean number of neutrons (curve) resulting from 
fission fragment of a particular mass number Af, with the 
associated dispersion indicated (bars), for 0.53 MeVn on 235 TJ 
(bottom) and 239 Pu (top). Data from Maslin [H, Nishio [3, 
and Tsuchiya [l5| are shown. [Color online.] 





V L 




V 


Clh 


n+ 239 Pu 


1.53 


1.43 


2.96 


-0.19 


n+ 235 U 


1.23 


1.23 


2.47 


-0.12 



TABLE I: The mean number of neutrons emitted from either 
the light fragment, Vl, the heavy fragment, Vh, or either 
fragment, V = Vl + vh, in fission events induced by thermal 
neutrons on U and Pu. Also shown is the correlation 
coefficient Clh = o(y L , vh)/\v(i>l)(t{vb)\- 
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Total number of neutrons v = v L + v H 



FIG. 5: Multiplicity distributions for neutrons emitted by the 
light (top), the heavy (middle), or either (bottom) fragment 
resulting from thermal-neutron induced fission of 235 U (cir- 
cles, blue) and 239 Pu (squares, red). [Color online.] 




Number of neutrons v 



FIG. 6: The mean total kinetic of energy the fission products 
together with the associated dispersions (bars), as a function 
of the neutron multiplicity in the event, for 0.53 MeVn on 
235 U (bottom) and 239 Pu (top). [Color online.] 



We note that the light fragment tends to emit more 
than its "fair share" of neutrons, a reflection of the fact 
that the excitation energy is not divided solely in propor- 
tion to mass. Futhermore, as the correlation coefficient 
Clh shows, there is a slight anticorrelation between vl 
and VH- This feature is presumably a result of the an- 
ticorrelation between the excitations of the two partner 
fragments caused by the thermal fluctuations of the heat 
partition at scission. 

Finally, Fig. [6] shows how the average total fragment 
kinetic energy of the fission products and their excitation 
depend on the number of evaporated neutrons v. The de- 
creasing character of the curves is easily understood since 
larger neutron multiplicities tend to arise from higher 
fragment excitations, which occurs in events with lower 
kinetic energies. 



V 


All 


1 


2 


3 


4 


5 


6 


7 




L 


2.30 


2.38 


2.30 


2.19 


2.02 








Pu 


H 


1.64 


1.70 


1.64 


1.58 


1.50 


1.34 


1.17 






L+H 


1.98 


2.10 


2.09 


2.01 


1.93 


1.82 


1.74 


1.68 




L 


2.18 


2.22 


2.17 


2.05 


1.85 








U 


H 


1.50 


1.56 


1.46 


1.39 


1.24 










L+H 


1.84 


1.85 


1.88 


1.84 


1.79 


1.73 


1.67 


1.55 



TABLE II: The mean kinetic energy e n (MeV) of the neutrons 
evaporated from the light fragment (L), the heavy fragment 
(H), or from either one (L+H), as a function of the respective 
multiplicity vl, vh, or v, in fission events induced by thermal 
neutrons on 235 U (bottom) and 239 Pu (top). 



C. Neutron energies 

We now turn to the kinetic energies of the evaporated 
neutrons. Figure [7] shows the fragment-mass dependence 
of the mean kinetic energy with respect to the frame of 
the emitting nucleus together with the associated disper- 
sion of the kinetic-energy distribution. 

The neutron spectra depend somewhat on the num- 
ber of neutrons emitted. This is summarized in Table |TT] 
which shows the mean kinetic energy of neutrons emit- 
ted from the light fragment, the heavy fragment, or from 




FIG. 7: The mean neutron energy e n (curves) together with 
its dispersion (bars) as a function of fragment mass Af for 
0.53 MeV n on 235 U (bottom) and 239 Pu (top). [Color online.] 
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FIG. 8: The mean kinetic energy and the associated disper- 
sions of all the neutrons emitted in fission events with a speci- 
fied total neutron multiplicity v, induced by thermal neutrons 
on 239 Pu and 235 TJ [solid curve, black dots), as well as the 
mean kinetic energy and the associated dispersions of all the 
neutrons emitted from the light fragment as a function of the 
corresponding multiplicity vl [dashed curve, green squares) 
and the mean kinetic energy and the associated dispersions 
of all the neutrons emitted from the heavy fragment as a func- 
tion of the corresponding multiplicity vh {dashed curve, red 
diamonds). [Color online.] 



either one, as a function of the respective neutron multi- 
plicities vl, vh, and v = vl + vh- 

The mean energies, as seen in the laboratory, as well as 
the associated dispersions, are displayed in Fig.[5]for the 
three neutron categories. In each case, there is an overall 
relatively modest decrease of the average neutron energy 
(and a corresponding narrowing of the distribution) as 
the neutron multiplicity is increased. This feature would 
be expected since the available energy must be shared 
among more neutrons. 

The full multiplicity-gated spectral shapes are shown 
in Figs. [9] (for U) and flOl (for Pu). It is apparent that 
the spectra become progressively softer at higher multi- 
plicities. This type of information is not provided by the 
standard models and is therefore novel. 



D. Neutron- neutron angular correlations 




Neutron kinetic energy e n (MeV) 



FIG. 9: The spectral shape of neutrons evaporated from the 
light (top), the heavy (middle), or either (bottom) fragment 
for specified values of the respective multiplicity vl, Vh, or v, 
in fission induced by thermal neutrons on 235 U. [Color online.] 




2 4 6 8 10 12 14 

Neutron kinetic energy e n (MeV) 
FIG. 10: Similar to Fig.® but for n+ 239 Pu. [Color online.] 



and 239 Pu. The analysis shown included only neutrons 
with kinetic energy above a threshold of 1 MeV. The 
results look qualitatively similar for other threshold en- 
ergies, with the angular modulation growing somewhat 
more pronounced as the threshold is raised (while the 
counting statistics is correspondingly reduced) . 

We see that the neutrons tend to be either forward or 
backward correlated. The backward correlation appears 
to be somewhat favored, as would be expected from the 
relatively small but negative value of the multiplicity cor- 
relation coefficient Clh shown in Table Q] 



The event-by-event calculation makes it straightfor- 
ward extract the angular correlation between two evap- 
orated neutrons, an observable that has long been of ex- 
perimental interest (see, for example, Refs. [H, [H, [2(| 
and references therein) but which cannot be addressed 
with the standard models of fission. 

Figure [TT] shows this quantity for the neutrons result- 
ing from fission induced by thermal neutrons on 235 U 



E. Neutron-photon correlations 

The final illustration is relevant for the correlation be- 
tween the neutron and photon multiplicities. Figure 1121 
shows the combined excitation left in the two product 
nuclei as a function of the total number of evaporated 
neutrons. When more neutrons are emitted the resid- 
ual product nuclei are less excited. This feature appears 



9 



CD 3 



J2 
CO 
_Q 

P 1 



-i 1 1 r- 



-i 1 r- 



-i 1 r- 



1 1 ■ I ■ i ' * 



60 90 120 

Relative n-n angle 9 12 
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Cn'y 


n+ 239 Pu 


2.97 


1.03 


5.67 


2.49 


-0.84 


-0.33 


n+ 235 U 


2.49 


0.96 


5.40 


2.40 


-0.79 


-0.34 



TABLE III: The mean neutron multiplicity V — (i>) and the 
associated dispersion <t„, the mean photon multiplicity 71 = 
(fi) and the associated dispersion cr M , and the neutron-photon 
multiplicity covariance ov M = (ufi) — vji together with the 
corresponding correlation coefficient C nl = ov^/fovo-p], for 
photons with transition energies above 200 keV. 



FIG. 11: The angular correlation between two evaporated 
neutrons for 0.53 MeVn on 235 U and 239 Pu, considering only 
neutrons with a kinetic energy above 1 MeV. [Color online.] 



CONCLUDING REMARKS 



to be reasonable since a larger-than-average number of 
neutrons is likely to have used up a larger-than-average 
portion of the total available excitation energy, thus leav- 
ing a less-than-average amount of residual excitation. 

Since the average number of photons emitted from a 
given product increases monotonically with excitation, 
the results in Fig. [12] provides a qualitative indication 
of the correlation between the number of neutrons evap- 
orated and the number of photons emitted during the 
further deexcitation of the product nuclei. Our simula- 
tions thus suggest that the two multiplicities are anticor- 
related: the more neutrons the fewer photons. 

This qualitative expectation is borne out by Tabic ITTT1 
which summarizes the result of including the actual pho- 
ton multiplicity /i into the analysis. The covariance be- 
tween v and /I is indeed negative and the corresponding 
correlation coefficient is about minus one third, suggest- 
ing a fairly significant degree of anticorrelation. 



t 1 1 1 1 1 1 r 




Total number of neutrons v 

FIG. 12: The mean total excitation of energy the two fis- 
sion products, together with the associated dispersions (bars), 
as a function of the neutron multiplicity in the event, for 
0.53 MeV n on 235 U (bottom) and 239 Pu (top). [Color online.] 



Over the last few years, experimental capabilities have 
improved dramatically while the practical applications of 
fission have broadened significantly. As a consequence, 
there has been an growing need for calculations of in- 
creasingly complex observables that are beyond the scope 
of the traditional models employed in the field. 

To meet this need, we have developed a new calcula- 
tional framework, FREYA, which can generate large sam- 
ples of individual fission events. From those it is then 
possible to extract any specific correlation observable of 
interest, without the need for further approximation. In 
developing FREYA, we have sought to make the numerics 
sufficiently fast to facilitate use of the code as a practical 
calculational tool. (Thus, on a MacBook laptop com- 
puter, it takes about 12 seconds to generate one million 
events.) 

Our early emphasis has been on creating a working 
code that can produce samples of reasonably realistic 
fission events and form a convenient basis for gradual 
refinements. (Its simple modular structure should fa- 
cilitate such further developments.) Consequently, the 
present version is still rather rough and cannot compete 
for quantative accuracy with established models, with- 
out suitable ad hoc parameter adjustments (see Ref. [|[). 
Even so, the model has already proven to be capable of 
making interesting predictions for correlations of interest 
in variety of contexts and we foresee an increased number 
of applications. 
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APPENDIX A: LIQUID-DROP MODEL 

For simplicity, we use here a the liquid-drop model plj 
for the macroscopic part of the nuclear binding energy. 
Accordingly, the surface and Coulomb energy of a spher- 
ical nucleus are given by 

2 



A 2 / 3 [l 



N - Z 
2A 



E° C ( A Z) = c 3 



A 1 / 3 



, (Al) 
(A2) 

17.9439 MeV, 



with the Lysekil parameter values: a 2 
k = 1.7826, c 3 = fe 2 /r = 0.7053 MeV [22|. 

The distortion energies of the prefragments at scission 
are based on the shape dependence of the surface and 
Coulomb energies of macroscopic prolate nuclei [f| 

(A3) 



E s (s) = E° c B s (e) 
E c (e) = E° c B s (e) 



E° s [l 

E° c [l-±s 2 ). (A4) 



45 



APPENDIX B: LEVEL DENSITIES 

The relationship between the nuclear excitation energy 
E* and the nuclear temperature T is generally somewhat 
complicated. For the tim e bein g, we simply use the famil- 
iar approximation T — e* /a, where e* = E* — 8E& e f is 
the statistical part of the excitation energy (the "heat" ) . 
For a given nucleus (Zj.Ai) the level-density parameter 
a,- is taken from Ref. 1231. 



ai(E* 



eo 



1 



SWj 
Ui 



(1 



Ui = E* 



(Bl) 

with eo = 7.25 MeV and 7 = 0.05. Here di = Ai/eo is the 
asymptotic level-density parameter whose parameter eo 
depends slightly on the specific value used for the damp- 
ing coefficient 7. The shell correction energies {<5Wi} and 
the pairing energies {A^} are those calculated by Koura 
et al. 10] for nuclei with 20 < Z t < 92. We note that 
a,{E* « A) « di{l + 5ff[l - \{E* - A,)]} is regular. 
Furthermore, a,i(E* = 0) ps Oj[l + r ySW l ] when jUi <C 1, 
which is most often the case. Finally, as E* is increased 
we have a,i(E*) — > a»[l + SWi/E*] — > a, = ^4i/e . 



APPENDIX C: SPECTRAL SAMPLING 

It is possible to devise a fast algorithm for sampling 
the spectral distribution (jTTj) for the evaporated neutron, 
dN/de oc ee~ e / T . It is based on the observation that the 
function xe~ x is a (normalized) Poisson distribution of 
order 2. Hence it can be expressed as the convolution of 
two (normalized) exponentials (each of which is a Poisson 
distribution of order 1), P2 = Pi * Pi with Pi(x) = e~ x , 

/>oc />oo 

xe~ x = / dxi I dx 2 8(xi + x 2 - x)e~ Xl e~ X2 . (CI) 
Jo Jo 

This is a special case of the general feature of Poisson 

P„ * Prt-n • 



We may therefore obtain a sampled value of the kinetic 
energy e as the sum of two energies, ei and e 2 , that have 
each been sampled from a usual exponential distribution 
oc e _£i / T . Since the sampling from an exponential distri- 
bution p(x) = e~ x is readily accomplished by sampling a 
random number rj that is uniformly distributed on the in- 
terval (0, 1] and then taking the negative of its logarithm, 
x = — In 77, the relative neutron kinetic energy is 

e n = ei+e 2 = -[\n m + \n m ]Tf^ , m G (0, 1] , (C2) 

where the two numbers rji have been sampled from (0,1]. 
Since both mean values and variances are additive un- 
der convolution and each exponential distribution yields 
(ej) = T and cr 2 (ej) = T 2 , the resulting relative kinetic 
energy e n has the mean value (e n ) = 2T™ ax and the vari- 
ance cr 2 (e n ) = 2(T^ ax ) 2 , for a fixed value of T™ ax . 

The energy spectrum of the post-evaporation photons 
can be sampled rapidly in an analogous manner, since 
the corresponding spectral shape, dN/de oc e 2 e~ £ / T/ , is 
(proportional to) a Poisson distribution of order 3, So 



-[In 771 + In 772 + ln7i 3 ]T/ , rj t £ (0, 1] 



(C3) 



It also follows that the mean value is (e 7 ) = 3T/ and the 
variance is c 2 (e 7 ) = 3T 2 , for a fixed value of Tf. 

APPENDIX D: LORENTZ BOOST 

We describe here the Lorentz boost required to express 
the motion of an ejectile and the corresponding daughter 
nucleus in the adopted reference frame. 

The boost velocity is that of the mother nucleus, Vj = 
Pi/Ei, where Pi is the momentum of the mother nucleus 
and Ei is its total energy, Ef = (M*) 2 + Pf. To perform 
the Lorentz boost, we first note that the component of 
the ejectile momentum parallel to the boost velocity is 
Pn = P n " v where v = Vi/Vi is the unit vector in the 
direction of V i. The component transverse to Vi is then 
pjj" = p n — pnV and this component is unaffected by the 
boost, = Pn- The parallel component of the ejectile 
momentum and its energy transform as follows, 

pi = 7(4 + E n Vi) , E a = 7 (E n + p n ■ Vi) , (Dl) 

where 7 2 = 1/(1 — V t 2 ) and E 2 — m 2 +p 2 . Thus the 
boosted ejectile momentum is 



Pn = [l E n 



7-1 
V 2 



Pn-Vi]Vi +p n , (D2) 



while the boosted value of the recoil momentum is ob- 
tained by reversing the direction of p, 



7-1 
V 2 



Pn -Vi]Vi -p n , (D3) 



distributions, P n 



with Ef being the total energy of the daughter, E 2 
(M* f ) 2 +p 2 . 



+m 
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